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ABSTRACT 


A  pre-stack  inversion  slgorithn  is  developed  for  sconstic  Kirchhoff, 
high-frequency,  con* on  offset  date.  Given  the  velocity  above  a  reflector, 
the  interface  is  located  and  an  angular ly-dependent  reflection  coefficient 
is  coaipnted  at  each  reflection  point.  A  quick  post-processing  step  then 
calculates  the  velocity  of  the  lover  nediun.  Latersl  velocity  variations  in 
the  second  layer  are  naturally  recovered  since  each  reflection  point 
provides  sn  independent  sieasure  of  the  reflection  coefficient.  The 
inversion  is  perfonaed  as  a  sapping  where  the  response  to  subsurface  test 
points  is  examined  by  an  integration  over  the  data.  If  a  test  point  is  on 
the  reflector,  the  reflection  coefficient  is  returned. 

Inversion  and  Migration  operators  both  utilize  an  integral  over  the 
data,  with  each  trace  in  the  suaaation  weighted  by  an  aaplitnde  and  a  phase 
tens.  Here,  knowledge  of  an  appropriate  inversion  phase  tens  is  gained  fro* 
a  Kirchhoff,  high-frequency,  forward  aodel.  To  deteraine  the  correct 
inversion  aaplitnde  term,  Kirchhoff  data  for  a  general  surface,  in  integral 
fora,  are  entered  into  the  inversion  operator.  The  resulting  integral  is 
evsluated  via  the  asyaptotic  aethod  of  4-diaensional  stationary  phase.  An 
aaplitnde  tera  is  then  chosen  so  thst  the  inversion  operator  produces  a 
singular  function  of  support  on  the  reflector,  weighted  by  the  reflection 
coefficient. 

The  Kirchhoff  offset  inversion  is  first  foraulsted  for  dets  acquired 
over  a  plane,  producing  a  3-D  reflectivity  asp.  Since  data  are  coaaonly 
collected  along  a  single  line,  a  2.5-D  specialization  is  slso  developed.  A 
aethod  for  deteraining  the  velocity  of  the  lower  aediua  froa  sn  angularly- 
dependent  reflection  coefficient  is  then  detailed  for  the  2. 3-D  case. 
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partial  derivative  matrix  of  phase  function  (22) 

inversion  amplitude  function  (14) 

compressional  velocity  (1) 

offset  variable  (B-10) 

free  apace  Green's  fnnction  (3) 

slope  of  2.S-D  inversion  output  interface  (58) 

source-receiver  midpoint  location  (13) 

unit  normal  vector  to  reflecting  surface  (7) 

cartesians  (2) 

test  point  location  (14) 

source  location  (2) 

receiver  location  (3) 

angular ly-dependent  reflection  coefficient  (6) 
source-reflector  distance  (16) 
reflector-receiver  distance  (16) 
source-test  point  distance  (13) 
test  point-receiver  distance  (13) 
tangent  vectors  to  reflecting  surface  (18) 

Kirehhoff  observed  field  ,  midpoint  coodinates  (15) 
observed  scattered  field  ,  midpoint  coordinates(13) 
backscattered  field  (44) 
scattered  field  doe  to  source  at  r+  (2) 
observed  scattered  field  (5) 

cartesians  for  backscatter  observation  point  (44) 
Dirak  delta  function  (3) 
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An  inversion  method  for  aeonstie  pre-stack  data  provides  a  means  for 
determining  subsurface  reflectivity  from  common  shot*  common  receiver,  or 
common  offset  time  sections.  In  exploration  for  hydrocarbons,  the  areas  of 
greatest  interest  are  those  where  the  geology  is  the  most  complicated.  As 
the  complexity  of  the  snbsnrface  increases,  the  validity  of  the  stacked 
section  as  an  interpretsble  depth  picture  decreases.  By  inverting  pre-stack 
data,  a  depth  reflectivity  profile  is  obtained  from  each  time  record. 

In  all  pre-stack  methods  the  reflection  coefficient  produced  is  for 
non-normal  incidence.  In  conaion  shot  or  common  receiver  gathers,  as  the 
shot-receiver  sepsration  increases,  so  does  the  angle  of  the  reflected 
energy.  Inversion  therefore  yields  the  reflection  coefficients  as  a 
function  of  angle,  and,  with  further  processing,  in  a  constant  density 
environment,  determines  velocities  below  the  reflector.  Analogously, 
processing  angulsrly-dependent  reflection  coefficient  data,  where  velocities 
are  known,  recovers  densities.  Common  offset  data  dictate  less  angular 
variation  of  the  reflection  coefficient  since  the  separation  between  source 
and  receiver  is  fixed  st  the  same  constant  distance  for  each  experiment. 

An  important  consideration  in  the  selection  of  an  inversion  data  set  is 
the  surface  acquisition  ares.  The  larger  the  data  collection  xone,  the 
greater  the  delineated  portion  of  the  subsurface  reflector.  In  this  paper, 
the  development  presented  is  for  the  inversion  of  common  offset  data,  for 
the  purpose  of  recovering  sub-reflector  velocities.  This  data  set  often 
contains  the  greatest  number  of  traces,  with  the  largest  areal  coverage. 

Inherent  in  the  definition  of  any  inverse  process  is  sn  assumption  of  a 
specific  forward  modeling  process.  Kirohhoff,  high-frequency,  non-zero 


offset  modeling  for  s  single  srbitrsry  surface  provides  the  basis  of  the 
inversion  developed  here.  Given  the  compressional  nave  speed  of  the  first 
medium,  the  vsve  field  effects  produced  by  this  Kirchhoff  date  are  inverted 
exactly  to  yield  the  interface  location  and  an  angnlarly-dependent 

reflection  coefficient.  Each  point  on  the  interface  has  a  unique  reflection 
coefficient  as  a  function  of  angle.  Since  the  offset  between  the  source  and 

receiver  is  constant  for  each  trace*  the  incident  angle  of  the  energy  to 

each  reflection  point  is  determined*  and  the  velocity  of  the  next  layer  is 
then  coaiputed.  Vith  each  reflection  point  independently  providing  the 
velocity  in  the  second  layer,  lateral  velocity  variations  beneath  the 

relector  are  naturally  recovered. 

As  with  migration,  inversion  is  performed  by  a  summation  of  traces, 
each  weighted  by  a  phase  and  an  amplitude  term.  The  inversion  phase  term  is 
of  opposite  sign  to  the  phase  used  in  the  Kirchhoff  forward  model.  The 
amplitude  term  is  an  unknown,  and  is  determined  by  inverting  lirchhoff  data, 
in  integral  form,  from  an  arbitrary  surface.  By  employing  the  asymptotic 
method  of  4-dimensional  stationary  phase,  the  inversion  operator,  acting  on 
high-frequency  lirchhoff  data,  is  evaluated.  The  unknown  amplitude  function 
is  then  selected  so  that  the  inversion  yields  a  singular  function  with 
support  on  the  reflector,  weighted  by  the  reflection  coefficient. 

The  Kirchhoff  offset  inversion  process  is  intuitively  similsr  to 
migration,  as  described  by  Schneider  (1978).  The  location  of  a  test  point 
is  input,  and  if  this  test  point  is  on  the  reflector,  the  value  of  the 
reflection  coefficient  is  returned.  For  each  test  point  an  integral  is 
performed  over  the  common  offset  traces,  with  each  trace  weighted 
geometrically.  The  weighted  traces  add  constructively  when  the  position  of 
the  test  point  coincides  with  the  reflecting  surface.  Since  data  are  not 


available  at  all  frequencies,  and  at  all  points  on  the  acquisition  plane,  a 
band-1 iai ted  and  apertnre-1 inited  singnlar  function  delineation  of  the 
reflector  is  prodaced,  as  deaonstrated  by  Bleistein,  Cohen,  and  Hagin 
( 198S) . 

The  inversion  is  exact  for  Kirchhoff,  high-frequency  data  froa  a  single 
interface.  Approxiaate  solutions  aay  be  obtained  in  aulti-layer  probleas  by 
a  single  pass  aethod  or  by  a  layer  stripping  technique.  The  accuracy  of 
single  pass  aulti-layer  inversion  is  dependent  on  the  input  velocity 
inforaation,  while  the  error  in  layer  stripping  inversion  is  a  function  of 
the  ability  to  downward  continue  the  wave  field,  preserving  aaplitudes. 
Neither  aethod  is  developed  here,  however,  for  coapleteness,  both  are 
described  in  the  context  of  the  derived  inversion  algoritha. 

The  inversion  algoritha  is  first  foranlated  for  3-diaensional  data  sets 
in  which  acquisition  is  over  a  plane.  For  this  data  base,  the  position  of  an 
arbitrary  3-D  interface  is  recovered,  along  with  an  angularly-dependent 
reflection  coefficient  at  each  reflection  point. 

The  algoritha  is  then  developed  for  the  2.S-D  case,  where  data  are 
acquired  only  along  a  line.  This  2.S-D  approximation  assuaes  invariance  of 
geologic  structure  in  the  off-line  direction,  thereby  introducing  a 
cylindrical  syaaetry  to  the  problea.  For  the  2.S-D  specialization,  a 
technique  is  presented  for  recovering  the  velocity  of  the  second  layer, 
given  a  reflection  coefficient  froa  the  inversion. 

Since  the  forward  aodel  is  of  such  fundaaental  importance  to  the 
inversion  development,  s  complete  Kirchhoff,  high-frequency  modeling 
derivation  proceeds  the  inversion  foraulation. 


KIBOnOFF,  lIGI-FtEXKJBNCl,  NON -Z  BIO  OFFSET  MO  DEL  INS 


For  the  pnrposes  of  inversion  end  siaulation  of  pre-steck  data,  a 
forward  aodeiing  procedure  is  required.  The  Kirchhoff  integral  aethod  is 
chosen  since  it  is  eabedded  in  wave  theory  and.  as  such,  produces  wave  field 
effects.  The  developaent  presented  here  is  suited  to  the  aodeiing  of  coaaon 


source,  coaaon  receiver,  or  coaaon  offset  gathers.  Of  particular  interest 
to  the  inversion  of  the  next  section  are  data  acquired  with  a  constant 
offset  between  the  source  and  receiver.  The  inversion  is  exact  with  respect 
to  this  Kirchhoff  forward  aodel  representation. 

A  high-frequency  assuaption  underlies  the  forward  aodeiing  and 
inversion  theory.  The  choice  of  a  suitably  high  frequency  is  a  function  of 
both  a  distance  paraaeter  of  the  problea,  and  the  velocity  of  the  aediua. 
By  assuaing  that  all  frequencies  of  the  data  are  greater  than  this  ainiaua 
high  frequency.  asyaptotic  evaluations  are  justified.  The  distance 
paraaeter,  r.  aay  be  the  ainiaua  depth  to  the  reflectors  of  interest  or  a 
"typical”  radius  of  curvature  of  a  reflecting  feature.  To  provide  a 
reasonable  approxiaation,  the  following  relationship  aust  be  satisfied: 


2wr/c  >>  1 


where  the  2,  appearing  in  equation  (1),  corresponds  to  the  2-way  travel  tiae 
of  the  forward  or  inverse  problea. 

For  exaaple,  if  the  depth  to  the  reflector  is  500  feet  in  s  median  with 
a  coapressional  wave  speed  of  10,000  ft/sec,  any  frequency  above  5  Bx  is 
considered  suitable  since  the  aaplitude  error  at  the  low  end  (5  Hz)  of  the 
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frequency  spectra  is  only  s  few  percent. 

The  snbseqnent  derivation  yields  the  high-frequency  lirchhoff 
representation  3>f  the  wave  field  reflected  from  an  arbitrary  surface. In  each 
experinent  the  source  is  offset  from  the  receiver.  For  modeling  purposes, 
this  experiment  is  repeated  along  the  surface  to  generate  a  time  section  of 
common  offset  traces. 

The  scattered  field  has  its  sources  on  the  reflector  and  is  thus 
governed  by  the  homogeneous  wave  equation  with  inhomogeneous  boundary 
conditions: 

V*Og(u,r,r+)  -  0  <2) 

The  two  spatial  variables  of  the  argument  of  0#  in  equation  (2)  indicate 
that  the  recorded  value  of  the  scattered  field,  0g,  at  any  point  r  is  a 
function  of  the  source  position  r+.  Since  the  scattered  field  is  recorded 
at  only  one  receiver  location,  r  ,  per  experiment,  a  sifting  operation  is 

performed  on  the  variable  r  of  equation  (2). 

For  the  purpose  of  sifting  under  a  volume  integral,  a  second  wave 

equation  is  introduced: 

V*g(«»  r,  r  )  +  w*/c  g(«,  r,  r  )  =  -  Mr  £  )  •  ^3) 

The  solution,  g(w, r, r"),  is  a  free  space  Green’s  function  which  describes  the 

propagation  of  a  point  source  from  the  location  r  ,  toward  any  point  r. 
Reciprocity  permits  r  and  r~  to  switch  places,  providing  another 


representation  of  equation  (3): 
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The  physical  situation  is  now  that  of  a  wave  traveling  from  any  point  source 
location,  r,  toward  a  receiver  location.  r~.  In  a  volume  integral  context, 
equation  (4)  describes  a  sifting  function  at  the  receiver  location,  and 
propagates  energy  from  a  reflector  to  the  receiver. 

The  wave  field  due  to  a  source  at  r+  and  a  receiver  at  r~,  U8(». r-,  r+) . 
is  obtained  in  the  following  manner:  multiply  equation  (2)  by  g(w, r“, r)  and 
equation  (4)  by  Df(w,  r,  r+)»  subtract  one  result  from  the  other,  integrate 
over  a  volume  of  space  bounded  by  the  reflector  which  contains  the  source 
and  receiver  points,  and  apply  Green's  second  theorem.  This  yields 


_  +  ff  f  +  ag(«.r,: 
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where  the  normal  vector  is  pointing  inward,  and  the  closed  surface  is 
composed  of  two  parts:  the  reflector  truncated  at  its  intersection  by  a 
large  hemisphere  having  its  base  on  the  reflector. 

The  scattered  field  is  given  exsctly  by  eqnstion  (5)  as  a  function  of 
frequency,  and  source  and  receiver  location.  To  compute  0#(w, r”. r+) ,  the 
scattered  field  and  its  normal  derivative  on  the  reflector  are  needed. 

An  approximate  solution  for  the  scattered  field  near  the  reflector,  due 


to  an  incident  point  source,  has  the  following  form: 


(6) 


exp[( iw/c) |r  -  r  |] 


0  (w,  r,  r  )  ~  8  . 

*  4n|r  -  r+| 


where  R  is  the  geometrical  optics,  angol arly-dependent  reflection 
coefficient.  A  derivation  of  this  high-frequency  reflection  coefficient  is 
presented  in  Appendix  A.  The  derivation  follows  that  of  Bleistein  (1984) 
and  is  included  for  the  reader's  convenience. 

The  high-frequency  approxiaation  of  the  noraal  derivative  of  the 
scattered  field  is  given  by 


80  (u,  r,  t+)  .  exp  [(iw/c)  |r  -  r+|l 

"  -  -  (iw/c)  (V|r  -  r  |*l)  R  - 
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The  free  space  Green's  function  solution  to  equation  (4)  is 


exp[( iw/c) |r  —  r 1 1 
g(w,  r  ,  r)  =  - 3 — : - - — 

4"l£_  -  £l 
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and  for  high  frequencies  its  noraal  derivative  is  approxiaated  by 


8g(  w,  r  ,r) 

5n 


~  (iw/c)(V|r  -  rj'ft) 
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Employing  the  Soanerfeld  radiation  condition,  the  integration  over  the 


hemisphere  is  neglected.  A  high  frequency  spproxiastion  for  the  scsttered 
field  is  now  obtained  f roa  equation  (S): 
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After  evaluating  the  gradient  terms,  the  scattered  field  due  to  an  offset 
source-receiver  experiment  is  written  as 
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where 
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Equation  (11)  provides  a  high-frequency  Kirchhoff  wave  field  response 
for  non-xero  offset  aodeling.  For  each  source-receiver  geoaetry,  a  surface 
integration  produces  a  single  trace.  Tiae  sections  are  then  constructed  by 
moving  the  source  and/or  receiver  in  the  desired  fashion. 


3-D  DiYBBSION  OPERATOR  DEVELOPMENT 


An  inversion  operator  it  developed  for  application  on  pre-stack  common 
offset  data,  prodncing  a  subsurface  nap  of  the  reflector,  along  with  its 
reflection  coefficients.  In  general,  a  reflector  ia  delineated  by  examining 
the  operator  response  of  many  different  subsurface  test  points.  For  each 
test  point  an  integral  is  performed  over  the  data  acquisition  surface,  with 
each  differential  contribution  of  the  integral  corresponding  to  a 
geometrically  weighted  trace.  When  a  test  point  is  on  the  reflector,  an 
sngolarly-dependent  reflection  coefficient  is  output.  The  velocity  below 
each  reflection  point  is  then  determined  by  a  quick  post-inversion 
processing  step.  Spatial  and  temporal  sampling  combined  with  finite 
acquisition  area  and  recording  time  dictate  a  band-limited  and  aperture- 
limited  singular  function  representation  of  the  interface.  The  peak 
amplitude  occurs  on  the  reflector  and  is  proportional  to  the  reflection 
strength. 

The  inversion  operator  that  is  derived  has  the  following  form: 


f[0B<«,m)]  -  f  |dm*A(  c,  «,  i >+.  R exp  [{-!•»/  c)  (■ ' 


The  frequency  domain  version  of  the  input  date  is  represented  by  Os(w.m), 
with  the  vector  m  parameterising  the  midpoint  associated  with  each 
experiment.  R*  +  is  the  distance  bvtween  the  source  and  the  test  point,  and 
R ’ ”  is  the  distance  from  the  test  point  to  the  receiver.  Ihe  singular 
function,  6(a),  acta  when  the  difference  between  the  test  point  position  and 
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the  true  reflection  position  equals  zero.  Pignre  1  illustrates  the 
psraaeters  of  the  problea,  with  variable  definitions  provided  in  Table  1. 

The  derivation  of  the  forward  aodeling  equation  provides  insight  as  to 
the  fora  of  the  inversion  operator.  In  the  forward  problea,  each  trace  is 
viewed  as  a  weighted  sub  of  iaage  sources  on  a  reflector.  Thus,  inversion 
of  forward  data  to  deteraine  the  interface  location  aust  involve  a  weighted 
sua  of  the  recorded  surface  traces.  Also,  the  idea  of  propagating  back  to 
sources,  in  aigration  or  inversion,  necessitates  an  inverse  phaee  tern. 
With  this  in  nind,  an  educated  guess  of  the  inversion  operator  is  aade: 

W[0#(a,a)]  “  Jjda2  jd«  ( -i«)B(a,  r,  r  \  c) exp[( -iw/c)  (R  ’+  +  R,_))U>(«»,b)  ,(14) 

where  B(a,  r,  r  ,  o)  is  an  unknown  function  containing  whatever  tens  were  left 
out  of  the  gness.  This  unknown  function  is  detemined  by  substituting  in 
known  analytic  Kirehhoff  data,  Uj^a.a),  for  a  single  interface  in  integral 
fora,  and  then  requiring  that  the  following  inversion  goal  be  realised 
asyaptotically : 

W[Ok(*,a)]  -  Ms)  R  .  (15) 

Upon  inserting  the  Xirohhoff  dsts  of  equation  (11)  into  the  inversion 
operator  defined  in  equation  (14),  the  following  relationship  is  obtained: 


W[Uk(*,«)l 


6(s)  I 


~ — 2  jjda2  |{d°2  jd“  [w2  B(g,  r,  r\e)  —  -  -+  - J  ? 

l'+  +  B  '“)]]J 


ezpl(iw/c) [(fi  +  B  )  -  (B 


(16) 


Note  that  the  differential  surface  element  of  equation  (11)  is 
functionally  represented  as 


dS 


l|T  dOj  do2 


(17) 


where  Oj,  o2  are  the  curvilinear  coordinates  that  paraaeterize  the 
reflecting  surface.  With  tj.  defined  as  the  tangent  vector  to  curves  of 
constant  o2,  and  t2  defined  as  the  tangent  vector  to  constant  oj  curves. 


I?  "  |*i  *  t2| 


(18) 


The  integral  in  equation  (Id)  is  evaluated  via  the  asyaptotic  method  of 
stationary  phase.  Since  the  forward  aodel  is  a  high-frequency  Kirchhoff 
representation,  there  is  no  further  assuaption  necessary  at  this  point.  In 
particular,  the  integrals  over  o  and  a  are  perforaed  by  4-diaensional 
stationary  phase,  with  the  integral  over  »  providing  the  singular  function 
of  arclength  along  a  curve  noraal  to  the  reflector.  The  unknown  function, 
B(a,  r,  r',  c)  is  then  chosen  such  that  the  relationahip  in  equation  (Id)  is 
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obtained: 


W[0,  ' w, ■)]  -  6(a) 


Given  an  integral  of  the  fora 


|  f(x)  exp[iag(x)]  dx  ,  x  ■  (xJf  Xj.  Xy  •••»  *  (19) 


containing  a  point  at  which  the  phaae  p(x)  is  stationary,  the  asyaptotic 
representation  for  large  a  is  given  by: 


K.)  -  pf]2  fU„> 


exp[iap(xQ)  +  i(sgn  A)n/4] 
'll  dot  A| 


The  point  of  atatioaazity,  xfl,  is  deterained  by  the  condition: 


▼*<In>  “  0 


The  aatrix  A,  at  the  stationary  point,  is  defined  as 


»Vso> 

j,»  - 1,2.3.  •••, 


sgn  A  =  2r  -  a 


where  r  is  the  nnaber  of  positive  eigenvalues  of  A.  a  is  the  order  of  the 


A\-v/vv-. 

*  O  O  O  a  4  <  * 


•  "  J-  „* 

^  ;• •  •  . 


/■v  V-ZvVv 


Vf 


i 

matrix,  and  m-r  ia  the  number  of  negative  eigenvr Idea.  The  difference  in 
the  nnaber  of  pocitive  and  negative  eigenvaloea  ia  therefore  equal  to  sgn  A. 
i  It  ia  shown  in  Appendix  C  that  no  eigenvalues  are  equal  to  zero  on  the 

reflector  for  the  Kirchhoff  inversion  problem.  Since  det  A  equals  the 
product  of  the  eigenvalues,  an  eigenvalue  of  zero  would  nullify  the  validity 
I  of  the  simple  stationary  point  asymptotic  integral  representation  of 

equation  (20). 

For  the  stationary  phase  evaluation  of  equation  (16),  the  phase  is 
taken  as 


*x) 


|r+  +  r]  -  [*'+  +  r ] 


(24) 


and  the  formal  large  parameter  is  w/c.  The  stationary  point  of  the  phase  is 
located  where  the  following  four  stationary  conditions  are  satisfied: 


do. 

x 


0 


i*l,  2 


(25) 


and 


(26) 


Details  concerning  the  derivation  of  each  partial  derivative  in  the 
stationary  conditions  and  in  the  matrix  A  are  provided  in  Appendix  B. 

The  conditions  of  stationarity  of  the  phase  are  as  follows: 
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"■"J I  ill  .1  f  jl « J .  J  J. 


VV.  'A  • 1 


t  Q'w  v  -  <  -  «  -  J 


-  0 


implies  J+  •  tj  -  8“  *  tj  .  1-1,2  , 


(27) 


U-.o 

8-l 


implies 


[  *i  ~  *‘i  -y]  -  hi  *  v  -  *i] 


[  i  ~  (,i  -  y]  -  [  (»i  *  di)  -  ill 


(28) 


Eqostion  (27)  is  s  font  of  Snell's  lav  which  reaffirms  that  important 
contributions  from  the  integral  come  from  specular  points. 

The  mstriz  A  is  expressed  as 


ill. 


a°i  a°j 


ill 


dm^  0Oj 


ill. 


doi  0,-j 


ih 


d\  aui 


(29) 


where  each  partial  derivative  is  a  2 -by- 2  matrix. 

In  the  asymptotic  evaluation  of  the  integral,  the  determinant  of  A  must 

be  calculated.  As  is  subsequently  demonstrated,  the  delta  function  of 

#+  + 

equation  (16)  acts  as  8  “  approaches  B~,  thereby  simplifying  the 
calculations.  This  is  analogous  to  the  test  point  approaching  the 
reflection  spot.  From  Appendix  B,  it  is  seen  that  when  -  R~, 
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i-1,2.  j-1,2 


(30) 


.v 


‘j 


0 


Therefore*  in  calculating  the  determinant  only  those  terms  which  do  not 
vanish  are  included.  yielding 


4|det  A  | 


det 


a!. 


8.,  8.J 


i*l»2,  j»l,2 


(31) 


It  can  he  demonstrated  from  equation  (B-29)  that 


J|det  A  |  -  z  |g 


(B*  +  B~)(B+2  +  B~2)(l  ~  8*  •  6~) 
(B+  B“)3  |6+  -  8'| 


(32) 


With  sgn  A  shown  to  be  equal  to  zero  in  Appendix  C,  the  asymptotic 
evaluation  of  equation  (16)  is  written  as 

W[U^(«, m)l  -  6(s)  B  - 

c2  B(m,  r,  r  , c)  B  [y  +  y  ]  (B+  B  )2  |8+  -  8  | 

4  z  (fi+  +  B~)  (B+2  +  fi"2)  (1  -  8+  •  8") 

fdw  exp{(ia>/c)[(B+  +  B")  -  (B>+  +  b'")])  .  (33) 


A  few  comments  are  now  required  to  describe  the  nature  of  the  function 


8(s).  In  a  1-dimensional  example*  6[p(t)]  is  given  ss 


equation  (34)  specifies  that  as  t  approaches  zero,  6[p(t)]  approaches  a 
spike,  weighted  by  l/|a|.  An  unweighted  delta  fnnction  is  therefore 
represented  as 


(35) 


In  the  3-diaensional  inversion  problen,  the  unknown  function, 
B(a,  r,  r  ,c),  is  obtained  such  that  the  inversion  operator  produces  a 
singular  function  scaled  only  by  the  reflection  coefficient.  The  weighting 
due  to  the  argument  of 


Mp)  -  S  [<e+  +  8")  -  (h'+  +  8'")] 


(36) 


■ust  therefore  be  accounted  for.  The  singular  function  of  equation  (36) 
acts  as  the  test  point  approaches  the  reflector  along  a  curve  noraal  to  the 
reflecting  surface.  The  weighting  factor  is  dp/dn,  and  an  unweighted 
singular  function,  6(s),  is  thus  deterained: 


with  the  ▼•riftble  s  denoting  ere  length. 

It  must  now  be  shown  that  there  is  only  one  point  at  which  the  delta 
function  acts,  where  R"  =  R~.  The  stationary  phase  conditions  from  the 
integration  over  the  data  acquisition  surface  are  rewritten  as 

sin  8*  +  sin  8^  =  sin  0j  +  +  sin  0j  (38) 

and 

sin  0j  +  sin  0^  =  sin  0^+  +  sin  0^+  ,  (39) 

with  the  equations  associated  with  constant  y  and  constant  x  planes, 

respectively.  Figure  2  illustrates  the  constant  y  case,  with  the  angles 

+  ,+ 

measured  counterclockwise  from  the  vertical,  and  the  vectors,  R~|  and  R-”j, 
defined  as  projections  onto  the  plane.  Note  that  if  the  test  point  is  to 
the  left  of  the  reflection  point,  sin  0*+j  <  sin  0+j  and  sin  0  <  sin  0“j, 

thereby  violating  the  stationary  condition  representation  of  equation  (38). 
Similarly,  the  test  point  cannot  be  to  the  right  of  the  reflection  point. 

It  must  therefore  be  inside  or  on  the  border  of  the  triangle  bounded  by  R+£ 
and  R-^.  An  analogous  argument  in  the  constant  x  plane  projection  indicates 
that  the  test  point  lies  inside  or  on  the  border  of  a  rectangular  cone  with 
an  apex  at  (x,y,z),  and  vertices  of  (m^-dj,  m2+d2. 0),  (mj-dj,  n^-dj,  0), 
(w,+d,,m,-d„,  0),  and  (b,^„b,^v  0) .  Finally,  equality  of  distances  within 


the  arguaent  of  the  delte  function  requires  thst  the  test  point  coincide 
with  the  reflection  point  for  the  delta  function  to  act. 


w 

In  order  to  solve  for  the  unknown  function.  B(a,  r,  r  ,  c), 
over  frequency  is  rewritten  in  the  following  fora: 

jdw  exp|( i«/ c)  [(R+  +  B  )  -  (8  +  +  R*  )]J  ■  2nch[(R+  +  R  )  - 

=  nc6(s)  |5*  -  8  | 
•+  •- 
(1  -  R  •  R  ) 


Noting  that 


,  +  .  -.  2(1"  •  R+  -  1) 

Iy  +  r  1  *■ 


.»+  ft-, 

i«  -  *  i 


the  value  of  the  unknown  function  is  now  obtained  as 


i*(b.  r,  r  ,  c) 


+  --+2  ,  »-2)(i  -  r+  .  8”) 


2z(R  +  R  ) (R  +  R _ 

,2  (R+  R-)2  |8+  -  8" | 


nc 


Inserting  this  function  into  equation  (14),  the  general 
pre-stack  inversion  foraula  is  produced: 


the  integral 


(R'+  +  r'-)] 

(40) 


(41) 


(42) 


3-D  Kirchhoff 
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“  H.  1  v*  *■  v  ■  *  ..  ■  .  *-.  -  .  •  .  ■ 


I[0#U.m)] 


~  6(s)  R  ~ 


*4  [LJ  L  „■  (,">  +  «">'»,,f2  *  »'~2)^(x  -  m ;  t- 

«cj ))  )  («•*  «-)i 

•  exp[(-iw/c)  (R,+  +  R  )]U  (u,a) 


(43) 


Since  the  singular  fnnction  acts  when  the  test  point  is  on  the 
+  |  + 

reflector*  R~~  is  replaced  by  B  The  angular ly-dependent  reflection 

coefficient  and  the  location  of  the  interface  are  therefore  determined  by 

eaploying  equation  (43)  on  coaaon  offset  data.  Within  a  depth  zone  of 

interest,  the  response  of  each  test  point  is  deterained.  If  the  input  point 

is  on  the  interface,  the  reflection  coefficient  is  returned.  Each 

differential  contribution  of  the  surface  integral  is  a  geometrically 

weighted  coaaon  offset  trace  in  the  frequency  doaain.  The  spatial  range  of 

integration  is  reduced  by  estiaating  the  aaziaua  dip  of  the  reflector. 

The  3-D  backscatter  inversion  result  is  obtained  as  a  special  case  when 

+  »  + 

the  offset,  d,  equals  zero  (R~  =  R  “) : 

W[Ot(w,£)]  =  jj^  jdw  o,  exp(-2io>R,+/c)  #0  (w,£)  =  Ms)  R  .  (44) 


This  checks  with  the  backscatter  result  of  Bleistein,  Cohen,  and  Hagin 
(198S).  In  particular,  after  an  integration  by  parts  is  perforaed  on 
equation  (42)  of  Bleistein,  et  al  (1983),  equation  (44),  above,  is  obtained. 
Ihe  necessary  integration  by  parts  is  described  in  Bleistein,  et  al  (1985) 
in  transforming  equation  (26)  to  equation  (27). 
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Extensions  of  the  method  to  the  more  realistic  case  of  a  mnl ti-layered 
earth  are  possible.  One  approach  is  that  of  layer  stripping  in  which  one 
"major”  reflector  is  inverted  for  with  each  pass  through  the  algorithm. 
Given  the  velocity  of  the  uppermost  layer,  the  initial  inversion  produces 
the  location  of  the  shallowest  reflector  and  the  velocity  of  the  next  layer. 
The  wave  field  is  then  downward  continued  from  the  acquisition  surface  to 
the  initial  reflector  with  a  procedure  such  as  that  demonstrated  by 
Berryhill  (1984).  This  initial  reflector  then  acts  as  the  acquisition 
surface  and  another  inversion  is  completed.  Problems  still  to  be  addressed 
in  this  procedure  include  that  of  amplitude  preservation  daring  the  downward 
continuation,  and  the  recognition  of  'suitable”  layers.  An  alternate 
approach  calls  for  a  single  pass  through  the  algorithm.  For  this  technique, 
the  velocity  employed  at  each  differential  integration  contribution  varies, 
corresponding  to  velocities  along  raypaths  connecting  each  test  point  to 
source  and  receiver  positions.  A  velocity  profile  is  therefore  required  as 
input  to  a  single  pass  multi-layer  Kirchhoff  inversion. 


2.5-D  INVEtSION  OPE1ATOS  SPECIALIZATION 


In  those  esses  where  dsta  collection  is  slong  s  single  line,  s  nr  it  l 
solntion  to  the  inverse  problem  is  not  possible  unless  off-line  reflector 
informstion  is  known.  An  sssnmption  commonly  made  is  that  the  direction  of 
data  acquisition  represents  the  direction  of  subsurface  geologic  variation. 
This  introduces  a  cylindrical  symmetry  to  the  problem,  indicating  that  all 
parallel  data  lines  would  produce  an  identical  time  section.  Under  this 
2.S-0  assumption,  the  inversion  double  integral  can  be  specialized  to  a 
single  spatial  integration  along  the  data  collection  line.  The  integral  in 
the  off-line,  invariant  direction  is  performed  by  the  asymptotic  method  of 
1-dimensional  stationary  phase. 

For  the  2.S-D  simplification,  the  data  are  taken  along  a  line  in  which 
the  coordinate  varies,  and  the  m2  coordinate  is  fixed.  The  reflector  is 
thus  assumed  to  be  a  function  of  x  only.  Snell's  law  for  this  cylindrical 
surface  dictates  that  the  y-value  of  all  specular  reflection  points  is 
identical.  Analogously,  in  the  inversion  operator,  the  important  part  of 
the  m2  integration  (the  stationary  point)  is  the  part  directly  over  the 
specular  reflection  line. 

With  data  collected  along  a  line  of  constant  m2,  *■>  integral  of  the 

following  form  is  considered: 


1(e)  =  jdm2  f (mj)  exp[iap(m2)] 


(45) 


The  asymptotic  representation  for  large  a  is 


1(a)  ~  f(y  )  - - —  exp{ia0(y')  +  ilsgn  (y' )\  n/4  (>gn  a)l  *(46) 

|o#  (y  )  |. 


where  y  ie  the  stationery  point. 

In  the  perticnler  case  of  equation  (43),  the  phase  is 


*<«,)  *  R*  +  +  R* 


and  the  stationary  point  satisfies 


0 


ie.  >2  s  y  •  The  second  derivative  of  the  phase  and  its  sign  at  the 
stationary  point  are  also  required  for  the  evaluation  of  the  integral: 


•*B  *  ly'  “  +1 


With  the  above  inforaation,  the  2.S-D  Kirchhoff  inversion  operator  is 


h' w  .  •  .  ■  .  *J*  *J»  V  */  V  ’«*  V  .* 


written  at 


I[0  (».»  )}  ~  iii  [da,  [d«  o>-^ - *■  *+  ■  j*  *  - -  ^2(1  -  8'+  •  8'  ) 

*  nc  J  1  J  (R  8  ) 


•  ezp[(-i»/c)  (R'  +  +  R*  )  +  (in/4)sgn(-u>/c)] 


2  rr _ 

-a  1  ~  j"" 

In  R 


1/2 


(51) 


Since 


u  «*  (sgn  «)  |u|  -  |u|ezp{( in/2)  -(in/2) (sgn  a)] 

■i |u|ezp[-( in/2) ( tgn  «)]  ,  (52) 

the  2.S-D  Kirchhoff  conaon  offset  inversion  operator  is  given  by 

W[U  («,■  )]  ~ 

s  1 

(1/2  _ 

f.  ,  ,1/2  (R<+  +  R '”)  <r'+  +  R'  )  J7I  8'  +  2". 

-■  r 1,1  — « — i(l  -  *  * 1 

•  ezp[(-iw/c)  (R,+  +  R*  )  -  (i3n/4)(sgn  »)]  D^w.Bj)  ,  (53) 

where  R>+  and  B  ~  are  evaluated  at  the  stationary  point: 


r  „  „,l/2 

R,+  =  [tz'  -  (.j  -  dj)]2  +  s'2]  ,  (54) 

r  ,  „l/2 

s'"  “  IKbj  +  dj)  -  z']2  +  s'2]  .  (55) 


For  each  test  point*  an  integration  is  perforaed  over  the  data  acquisition 


line 


If  the  teet  point  is  on  the  reflector,  the  angular ly-dependent 


reflection  coefficient  is  output. 

i  t 

The  inversion  produces  s  reflectivity-depth  section,  where  each  x  ,x 
interface  coordinate  has  a  unique  reflection  coefficient  associated  with  its 
unique  incident  energy  angle.  Computation  of  the  lower  layer  velocity  at 
any  reflection  point  necessitates  the  measurement  of  the  interface  slope, 
h'(x),  from  the  inversion  section. 

With  the  upper  layer  velocity,  c,  known,  and  denaities  assumed 
constant.  Appendix  A  provides  a  relationship  for  the  lower  layer  velocity, 

cl: 

r  ,  « 

|  2  2  .  -21 

C1  “  lYt  "  Yi  +  *  J  *  (56) 

where 


Ti 


i+ 


c 


( A-26) 


and 


Tt 


(1  -  8) 
7TT-RT 


(57) 


The  normal,  ft,  is  given  by 


h  *  (x)  f  -  £ 
[.  * 


(58) 


leading  to 


napped  by  deteraining  the  lower  Bedim  velocity  f roa  a  range  of  reflection 
pointe  along  the  interface. 

To  deteraine  velocities  froa  a  3-D  inversion  output.  slopes  are 
aeasnred  in  two  orthogonal  directions.  Two  relationships  siailar  to 
equation  (61)  are  produced,  in  two  unknowns,  Bj  and  Once  a^  and  B2  are 

deterained.  c.  is  obtained,  ea ploying  a  3-coaponent  noraal  vector. 


A  pre-stack  inversion  operator  has  been  developed  for  common 
offset  data.  Given  the  velocity  above  the  reflector,  the  interface 
location  and  its  reflection  coefficients  are  determined.  A  data  set 
acquired  over  a  plane  yields  a  3-D  reflectivity  map  of  the  interface. 
When  data  are  available  along  a  line,  it  is  assumed  that  the  data 
collection  is  in  the  direction  of  geologic  variation.  A  2.S-D 

cylindrioally  symmetric  reflectivity  mapping  is  then  obtained.  Given 
reflection  coefficients,  a  quick  post-processing  step  determines  sub¬ 
reflector  velocities  after  measuring  interface  alopea  off  of  the 
inversion  output  depth  sections. 

The  inversion  is  exact  for  Kirchhoff,  high  frequency,  common 
offset  data  from  experiments  over  a  single  reflector.  If  the  reflector 
is  S00  feet  deep  in  a  medium  with  a  compressional  wave  speed  of  10.000 
ft/sec,  any  frequency  information  greater  than  5  Hx  is  suitable  as 
input. 

There  are  several  extensions  to  the  single  interface  inveraion 

method  that  could  address  the  multi-layer  problem.  The  following 
suggested  idess,  however,  neglect  multiples.  One  approaoh  is  to  assume 
that  each  layer  is  independent  of  all  other  layers.  In  this  esse  the 
inversion  errors  are  sensitive  to  the  choice  of  e  single  "  above 
reflector”  velocity.  Another  method  is  to  provide  multi-layer  velocity 
information,  and  then  employ  Snell's  law.  Travel  times  along  raypaths 


connecting  each  test  point  with  nil  source  and  receiver  positions  are 
then  computed.  This,  however,  requires  knowledge  of  subsurface 
velocities.  A  final  idea  is  that  of  layer  stripping.  The  initial 
inversion  yields  the  location  of  the  first  reflector.  A  wave  equation 
downward  continuation  of  the  data  is  then  performed  to  this  reflector. 
The  subsequent  inversion  velocity  is  determined  from  the  reflection 
coefficient  and  an  inversion  to  the  next  reflector  is  completed.  The 
layer  stripping  method,  however,  assumes  that  the  data  are  of  such 
quality  that  layers  are  discernible  across  the  section.  Also,  a 
downward  continuation  procedure  is  required  for  pre-stack  data  which 
preserves  amplitudes. 
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APPENDIX  A:  AMBULAILT-DBPBNDfiNT  REFLECTION  COEFFICIENT 


The  angularly-dependent  reflection  coefficient,  8,  ie  determined  nnder 
the  assumption  that  the  surface  appears  locally  planar  to  incoming  and 
outgoing  energy.  This  is  a  high-f reqnency  approximation  which  permits  title 
calculation  of  8  via  plane  wave  analysia,  nnder  the  conditions  of  continuity 
of  the  field  and  its  normal  derivative  across  the  interface. 

The  frequency  domain  representation  of  the  incident,  reflected,  and 
transmitted  field  is  given  by: 

Oj  ~  A^expduv^)  ,  ( A— 1 ) 

0^  ~  A^expdwc^)  ,  (A-2) 

Ut  ~  Atexpdut  )  .  (A-3) 

Each  amplitude  term.  A,  is  a  function  of  space  and  frequency,  and  t  is  the 

eikonal  function,  equalling  travel  time  along  the  raypath. 

The  total  field  above  the  reflector  is  given  as 

0^  ■  A^expdwtj)  +  A^expdwc^)  ,  (A-4) 

and  the  total  field  below  as 

0  »  A  expdmc  )  .  ( A— 5 ) 

DC  t 

Demanding  continuity  of  the  field  on  the  reflecting  surface  requires  that 

A^expdtoTj)  +  A^expdwt^)  =  A^xpdw^)  ,  (A-6) 


which  is  possible  only  if  the  phases  match  on  the  surface.  Therefore  the 
first  continuity  condition  yields 


Continuity  of  the  norael  derivative  of  the  field  across  the  surface 

states  that 

(VAt  +  iwAjVr^-8  +  <YAc  +  iuArVTr)*8  =  (VAt  +  iwAtVTt)*8 

.  ( A— 8) 

The  high  frequency  assumption  coupled  with  the  relatively  small 

of  amplitude  yields: 

variations 

Vi 4  V. '  Vi  • 

where 

( A-9) 

Tj  - 

( A-10) 

Tf  *  TV*  • 

( A— 1 1 ) 

rt  *  ▼«t** 

( A-12) 

Employing  equations  (A-7)  and  (A-9),  At  and  At  are  solved  for  in 

and  the  eikonal  gradients: 

terms  of  Aj 

Yi  “  Tt 

A  =  — - —  A, 

r  -rt  ♦  r,  i 

( A-13 ) 

Ti  "  Tr 

*«  "  -rt  *  r,  *.  • 

( A-14) 

It  is  now  necessary  to  solve  for  yf  and  yt  in  terms  of  y^. 

above,  the  phases  must  match  on  the  boundary. 

As  stated 

5 


l  r  t 

Taking  the  tangential  derivative  of  equation  (A-1S)  produces 

Vt^tj  -  VTf*tj  -  Vt^tj  ,j  -  1,2  ,  ( A-16) 


where  tj  and  tj  represent  two  unit  vectors  which  paraaeterize  the  surface. 
Since  the  tangential  coaponents  of  the  eikonal  function  gradients  aatch  on 
the  interface,  the  total  gradients  on  the  interface  are: 


Vti  -  T  +  Yil 
Vtr  -  T  +  yrft 
Vrt  -  T  +  yt* 

where 


Substituting  equation  (A-l),  (A-2) 
equation  yields  the  eikonal  equation 


( A-17) 
(A-l  8) 
(A-l 9) 


tj)  tj  .  ( A-20) 

or  (A-3)  into  the  hoaogeneous  wave 


|Vt|  -  1/c 


( A-21) 


Therefore,  dotting  equations  (A-17)  and  (A-18)  with  Vt^  and  Vtf, 
respectively,  gives 


l^tj2  =  |T|2  +  Yi2  =  c-2  (A-22) 

|VtJ2  *  |T|2  +  Tt2  -  c"2  .  ( A-23) 

Since  the  incident  and  reflected  energy  are  of  opposite  sense  with  respect 
to  the  surface  normal  direction,  the  following  relationship  is  obtained: 

Yj  -  "  Tr  •  ( A-24) 

Analogous  operations  to  equations  (A-18)  and  (A-19)  provide 

Tt  =  ~  +  Y2  ,  ( A— 2 5 ) 


where  the  velocity  in  the  upper  medium  is  c,  and  the  lower  medium  velocity 
is  Cj .  For  the  current  consideration  of  a  two  layer,  piecewise  constant 
velocity  medium 


Yj  =  (8*8  )/c 


( A-26) 


where  R+  is 


the  unit  vector  along  the  raypath  from  the  source  point  to  the 


reflector. 


t 


r4  +  rt 


i 


Taking  the  ratio  of  the  reflected  to  the  incident  amplitudes  produces 
reflection  coefficient: 


B 


Similarly,  the  transmission  coefficient  is 


\ 

\ 

) 

C  A— 29) 

i 

* 

i 

i 


Ti  +  rt 


<  A-30) 


t 

I 


( 


» 


i 

I 
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APPWDII  B:  STATIONARY  PHASE  CALCULATIONS 


Coaputation  of 


The  derivative  of  the  phaae  function  with  reapect  to  the  paraaietera 
defining  the  reflector  ia  given  by 


d 


[b+  +  b~]  -  [b'+  +  b'“]  J 


i-1.2 


(B-l) 


lx  lx 

Since  B  and  fi  are  functions  of  an  independent  teat  point  and  therefore 
unrelated  to  the  true  reflector,  the  partial  derivative  reducea  to 


%  '  5§7  [,+  *  *’]  i-l.: 


(B-2) 


Expanding  the  first  term  results  in 


dB+  .  ~+ 

$o^"  dx  IcT 


dfi+  dx  ^  dfi+  dy  *  dB  dz 

7T-JT  *  JTTT  TT  TT 


i-1.2  , 


(B-3) 


which  is  equivalent  to 


dB+ 


VB+  •  t.  =  8+  •  t.  .i-1.2 


-i 


ii 


(B— 4) 


Analogously, 
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Therefore,  the  partial  derivative  of  the  phaae  with  reapeet  to  the 


variables  which  paraaeterize  the  reflector  ia 


*i 


,  1=1,2 


(B-6) 


The  firat  condition  of  stationarity  of  the  phaae  ia  expressed  as 


»  i-1.2 


(B— 7 ) 


or 


fi+  •  t£  -  1“  •  tt  ,  i=l. 2  .  (B-8) 

This  relationship  states  that  when  evalnating  the  integral  of  equation  (17), 
the  iaiportant  contributions  cone  froa  the  apecular  points  associated  with 
Snell's  law. 


Coaputation  of 


The  derivative  of  the  phase  function  with  respect  to  the  aidpoint 


variables  of  the  data  acquiaition  aurface  is  given  by 


i-1.2 


(B-9) 


^  [  \r  -  K  *  •■] 


Examination 

of  the 

firat 

two 

tenia  providea 

insight  to 

the 

entire 

derivative: 

3R+ 

lx  - 

("l  " 

VI 

3R+ 

“ 

ty  -  <«2  - 

d2)J 

R+ 

«♦ 

3R~ 

[(■j  + 

V  " 

i] 

3R~ 

"*!  +  V 

-  y] 

(B— 1 0) 

9-l 

R' 

- 

0-2  " 

R_ 

““  “  n 

In  •  snccint  notation,  the  entire  partial  derivative  ia  written  aa 


l£r 


[Xj  -  <«j  -  ) ]  [(■i  +  d^  -  x^ 


i  B  R 

lx’  -  (Wj  -  dj)J  +  dj)  -  xj  ] 


, »  + 


,  i-1.2  .  (B-H ) 


where  Xj  -  x,  and  x*  “  y.  The  aecond  condition  of  atationarity  of  the  phaae 


ii.li  ill 1 1  .115.111. il « 


i-1.2 


(B-12) 


or 


I*  .  t  -  I'  •  t  -  I".l-I'  !t 

r .  f  - 1  •  j  -  i~ .  j  -  r-  •  j 


(B-1S) 

(1-14) 


Computation  of 


do^do 


j 


The  second  partial  derivative  of  the  phase  with  respect  to  reflector 
parameters  is 


a.,  a.j 


,  i-1.2  a  j-1.2 


(B-15) 


Since  t ^  varies  as  a  function  of  o^ 


where 


1  i-J 
0  i*j 


(B-16) 


(B-17) 


Therefore; 


Separating  the  renal ning  decivativea  into  two  terns  yields 


a8+ 


a 


ii 


a 


».  i.  •  ?  «ij  •  *i>  -  7  <**  •  *i»<r  •  it> 

J  K  K 


and 


3R 

3o" 


af]  i  , 

-  •  t.  -  y*;  -  •  t.  - - ft.  •  t.)  +  —  (S  •  t.)(«  •  t4>  .(B- 

j  3j  ~l  i~  B  “j  1 


20) 


Conbining  all  tema.  the  aeeond  partial  derivative  is  obtained: 


-  7  ta  •  *»* " <!*  •  v<s+  •  ‘*’1 

♦  [<£,  •  £,>  -  (i  •  £j><»-  •  £,)] 

* (,t  ‘  r>  •  sij  *tj  • 


(B-21) 


Conpntation  of 


±il 


The  next  second  partial  derivative  is  taken  with  respect  to  both  the 


acqnisition  surface  Midpoint  variables,  and  the  reflector  paraaeters. 


,  >2>  .  f(t*  -  I-)  .  ,  1 

*'71=7  I'  £ij 


i-1,2 


j-1.2 


(B-22) 


Since  the  two  sorfaeet  ere  independent  of  one  another  the  partial  derivative 
ia  redneed  to 


a"oj  a«j  "  fclj***  "  8_)]  *  ^i  1‘1*2  *  J"1*2 


(B-23) 


Thia  firat  ten  ia  evalnated  aa  follows: 


ar 


dmi  -i 


t.  - 


J-1.2 


t_  !L .  i_  t+l .  t 

U*  *"j  t+*«,  J  11 

i_  [., .  b<  :.p :  V1  i*  [  •  st  w 


(B-24) 


(B-25) 


where  1  is  the  on  it  vector  in  the  x-direction  and  !  ia  the  unit  vector  in 
the  y-direction.  The  second  ten  of  the  partial  derivative  ia  aiailarl7 
obtained: 


£ 

R 


lii 


R 


*i 


J-1.2 


Combining  these  two  tens  and  noting  that 


(B-26) 
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. -  .  "  \  -  .  *.  l".  ^ ’V*  ‘  ■  a  ■«  *  .  * 


ft: 


[<«  +  d  )  -  X  ] 

- * - 1 - i-  -  8  •  J 


(B-28) 


the  final  fore  of  the  second  partial  derivative  is 


>\  . 

[  t  *  rF’  *  ' ,,l+  *  r1’’ ,)8’  ] ' £i  «■*.  J->- 


2  .  (B-29) 


A 

Computation  of  i- 

dmi  d.j 

The  final  second  partial  derivative  is 


a-i  dmj 


l_  f.  [*1  -  (>1  -  Ji>l  +  [(,i  *  V  ~  *l] 

*J  l  I*  «' 


+  K  -  (,i :  v] 


[(«t  +  at)  -  «;  ] 


1  1-1,2.  i-1.2.  (B-30) 


Each  term  is  obtained  in  the  saae  Banner*  the  first  is  illustrated: 


AFMMDIX  C:  S6N  A  CALCULATION 


The  value  of  det  A  end  tgn  A  ere  required  et  the  point  of  action  of  the 
singular  function.  At  this  point  it  is  seen  ,  fron  equations  (29)  and  (30), 
that  det  A  is  non-negative,  regardless  of  the  shape  of  the  reflector. 
Furtheraore,  since  8*  and  8”  are  of  opposite  sense,  as  depicted  in  Figure  1, 
equation  (32)  states  that  det  A  is  never  equal  to  zero  on  the  reflector. 
Since  det  A  equals  the  product  of  the  eigenvalues,  this  product,  X2X2X3X4, 
is  always  positive  on  the  reflector.  The  sign  of  the  eigenvalue  product 
reaains  unchanged  if  the  sign  of  any  two,  or  of  all  four,  eigenvalues 
changes.  However,  to  change  sign,  an  eigenvalue  aust  pass  through  zero,  at 
which  point  det  A  equals  zero.  Therefore,  since  det  A,  on  the  reflector,  is 
never  equal  to  zero,  the  sign  of  each  eigenvalue  is  constant  on  the 
reflector,  and  is  not  a  function  of  its  shape.  Vith  the  shape  of  the 
reflector  not  influencing  the  eigenvalue  signs,  the  siaplest  interface,  a 
flat  plane,  is  employed  in  calculating  sgn  A. 

For  a  flat  plane,  let  tj  =  f,  £2  *  J»  Oj  “  z,  and  03  “  y.  Without 
loss  of  generality,  data  acquisition  is  in  the  z-direction,  vith  d2  *  0 
and  dj  ■  d.  The  stationary  conditions  specify  Bj  <■  z  and  B2  ~  y*  yielding 


8+ 

8- 


1  -  R 

P 

dt  +  zk 

- S - 

p 

df  -  zk 


[d2  ♦  Z2] 


1/2 


(C-l) 
( C-2) 

( C-3) 


The  aatriz  A  is  now  written  as 


To  siaplify  the  algebra  involved  in  aolving  the  characteristic  equation  of 
aatrix  A,  let  the  following  definitions  apply: 


(C-5) 

(C-«) 

(C-7) 


The  eigenvalues  are  obtained  by  solving  the  algebraic  equation  of  degree  4, 

I 


det  (A  -  XI)  =  0 


( C-8) 


or  when  expanded: 

X4  -  (a  +  c)X3  +  (ac  -  c2  -  b2)X2  +  (ac2+  cb2)X  +  <cb)2  -  0  .  (C-9) 


Sgn  A  is  deterained  froa  equation  (C-9)  without  explicitly  finding  the 
eigenvalues.  The  product  of  the  eigenvalues  is  (ob)*.  Since  this  product 


is  always  positive,  no  eigenvalues  are  equal  to  zero.  Furthermore,  this 
restricts  the  signs  of  the  eigenvalues  to  three  cases:  all  are  positive,  or 
all  are  negative,  or  there  are  two  of  each  sign.  The  sum  of  the  eigenvalues 
is  (a+c),  which  is  also  positive.  Therefore,  not  all  of  the  eigenvalues  are 
negative.  Finally, 

-(.c2  *  ct>2)  -  Jj  s 

with  (ac2  +  cb2)  always  positive,  indicates  that  some  of  the  eigenvaluea  are 
negative.  Thus,  two  are  positive,  two  are  negative,  and  sgn  A  equals  zero. 

It  is  possible  for  sgn  A  to  be  nonzero  far  from  the  reflector.  In  the 
2-diaensional  stationary  phase  calculations  perforaed  in  Eirchhoff  forward 
■odeling,  Cohen  and  Bleistein  (1984)  demonstrate  that  passing  through  buried 
foci  changes  sgn  A.  The  resulting  distributions  have  support  on  the 
reflector,  with  negligible  amplitudes  far  froa  the  reflector.  Analogously, 
foci  associated  with  the  2  aurfacea  of  the  inversion  problem  aay  produce  a 
nonzero  sgn  A.  However,  the  large  paraaeter  of  the  asymptotic  analysis 
requires  large  radii  of  curvature,  placing  foci  far  from  the  reflector. 
Therefore  contributions  of  these  nonzero  sgn  A  distributions  are  negligible. 


r+i  k=j  +i  1  1 


(C-10) 


TABLE  1 


Integration  Parameter  Definitiona 

source  location  :  r+  =  (£+,  0)  =  (?+,  n+,  0)  »  $+f  +  tj+J 
receiver  location:  r”  =  (£",  0)  =  0)  ■  {"t  +  r|~  J 

sonrce-receiver  midpoint:  m  =  (mj,a2*0) 

offset  variable:  d  «=  (dj.djiO) 

source  location  in  midpoint  variables:  r +  *  (rn^-dj. mj-d2» 0) 
receiver  location  in  midpoint  variables:  r~  *=  (mj+dj, m2+d2»  0) 

reflection  point  location:  r  =  (x,y,z) 

sonrce-reflector  vector:  B+  *  r  -  t+  *  <x-$+, y-n+» a) 
reflector-receiver  vector:  H_  «  r-  -  r  *  ( $_-x,  n~-y,  -z) 

teat  point:  t’  *  (x*,y  .s’) 

source-test  point  vector:  B  «  r  -  r  *  (x  -{,y  -r\ §  z  ) 
test  point-receiver  vector:  R  =  r  -  r  -  x  ,  n  -y  , -i 
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